Numerical Study of the Three-Dimensional Gauge Glass Model 
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We investigate numerically the finite-size scaling properties of the domain wall energies in the 
three-dimensional gauge glass model. From the analysis of results obtained for systems of linear 
sizes 3 < 1/ < 8 we conclude that the stiffness exponent of the model is positive. This implies the 
existence of a stable ordered phase at low but finite temperatures. 



PACS numbers: 64.60. Cn, 75.10.Nr, 05.70.Jk 

It has been suggested that in type-II superconductors 
at low temperatures defects may pin the flux lines at ran- 
dom positions thus destroying the Abrikosov vortex lat- 
tice. This leadS|-tfl a new type of superconducting state, 
the vortex glasxM, in which the phase of the supercon- 
ducting order parameter is random in space but frozen 
in time, much in the same way magnetic moments are 
frozen in the low-temperature phase of spin glasses. The 
simplest system expected to have an ordered phase anal- 
ogous to the vortex glass is the gauge glass model, origi- 
nally introduced to describe disordered arjays of Joseph- 
son junctions in an external magnetic fieldl3'cl. This model 
is defined by the Hamiltonian 
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where 9i is the phase of the order parameter at the i- 
th site of a simple cubic lattice and the sum runs over 
all pairs of neighboring sites. The energy scale is set by 
the coupling constant J and the lattice spacing is iden- 
tified with the typical distance between vorticesa. The 
phase shifts Aij = (27r/$o) // A-dl where A is the vector 
potential of the applied magnetic field and <I>o the flux 
quantum. The effects of the disorder in the positions of 
the vortices are incorporated by taking the phase shifts 
as independent quenched random variables. The situa- 
tion that interests us, where the disorder and the external 
field are large, may be modeled by taking Aij uniformly 
distributed in the interval [0, 27r]. The three-dimensional 
gauge glass model has been exteHsively studied numeri- 
cally by Monte Carlo sinmlationljlS and finite-size scal- 
ing of defect wall energiesuclM the important issue being 
whether a thermodynamically ordered phase can exist at 
finite temperature in this system. pjAlthough the results 
of the earlier Monte-Carlo studiesOu of the model were 
consistent with the existence of a low-temperature vortex 
glass phase, they could not rule out a zero-temperature 
transition since only small systems could be brought to 
equilibrium below T ^ 0.6J. Stronger evidence in favor 
of a finite-temperature transition has been obtained in 
recent simulations based on the vortex representation of 
the problemQ in which it was found that the transition 
temperature may be as high as Tc = 0.93 J. On the other 
hand, the domain-wall renormalization-group (DWRG) 



studies performed so faiil0 were inconclusive, the sizes 
of the systems studied being too small and the statisti- 
cal error too large to decide unambiguously whether the 
lower critical dimension of the model is above or below 
three. In this paper we reexamine this problem by means 
of a DWRG study of model ([^) with an algorithm that we 
have recently proposed and-applied to the XY spin-glass 
model in three dimensionsEll. This algorithm allows us to 
study lattices substantially bigger than with conventional 
methods as well asijpJmprove upon the statistics. In the 
defect wall methodtaa the energy cost W of introducing 
a domain wall in the system is studied as sufpjction its 
linear size L. In the scaling regime one findaij'EJ W ^ LP 
where the stiffness exponent d may be positive or nega- 
tive depending on whether the system is above or below 
its lower critical dimension . From the results obtained 
for our five largest sizes (4 < L < 8) we find the value 
^'gg = 0.077 ± 0.011 for the gauge glass model. This re- 
sult implies that a stable ordered phase exists at low but 
finite temperatures. 

To determine the domain-wall energy one computes 
the differences Ai? — Ep — Ea between the ground-state 
energies corresponding to periodic (P) and anti-periodic 
(AP) boundary conditions along some direction for an 
ensemble of systems of size Ng = L^. The boundary 
conditions along the two remaining directions are kept 
fixed. For sufficiently large systems the distribution of 
energy differences diffierences V{AE,L) is expected to 
have the scaling forn£3 



ViAE,L) = L~'^V{AEL-'^). 
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The width of the distribution, W{L) = {AE^)\^ 
is interpreted as tlie [effective coupling constant between 
blocks of Ng sitesErEJ. If 6* > 0, the rigidity of a block 
diverges with its size, which indicates that the system has 
long-range order. If the stiffness exponent is negative, the 
correlati(jai| length diverges at T = with ^ ~ T~'^ and 

The ground state of the gauge glass model is given by 
the absolute minimum of ([^) subject to the appropriate 
boundary conditions. In the presence of disorder the ex- 
tremal conditions dH/dOi = Vi have in general a very 
large number of solutions whose presence greatly compli- 
cates the task of searching that with the lowest energy. In 
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the spin-quench algorithmlli(SQA) usually employed to 
solve this type of problem, long sequences of metastable 
states are randomly generated among which one will find 
the ground-state provided the number of trials is suffi- 
ciently large. Since the number of metastable states o£a 
frustrated system increases exponentially with its sizdU, 
so does the number of trials required. This limits the 
maximum size of the systems that can be studied us- 
ing this method in practice. We have recently proposed 
a far jaore efficient algorithm for the search of ground 
stateal^. It is based upon the morphological characteris- 
tics of the low-lying states of frustrated XY models as m- 
vealed by detailed examination of numerous examplesEj. 
For a given realization of the disorder in Eq.|l|, the low- 
energy configurations are characterized by the existence 
of regions where the order parameter varies smoothly (do- 
mains), and others where the spatial distribution of the 
phases looks pretty much random. The former exist in 
parts of the sample where frustratian is low, the latter 
where it is high. As it turns outll3, the position, size 
and shape of the domains are mostly determined by the 
realization of the disorder and are essentially the same 
for all the low-energy states. Aside from smooth dis- 
tortions of the order parameter, the essential differences 
between any two such states are almost rigid rotations of 
the individual domains, accompanied by large amplitude 
rearrangements of the phases in the frustrated regions 
between them. Stationary states in which the domain 
structure is disrupted do exist, but their energy is much 
higher. In our method, sequences of low-energy configu- 
rations are generated recursively in such a way that the 
domain structure is preserved at each step. The result is a 
reduction of the probability of appearance of high-energy 
configurations in the sequence and a corresponding en- 
hancement of that of finding the ground-state or states 
lying nearby in energy. The procedure is as foUowgiil. 
The first state in the sequence, {fl*-"-'}, is obtained by 
a conjugate-gradient minimization (CGM) of the energy 
(|l|) starting from a random distribution of phases. New 
states are generated by iterating the following steps, i) 
Sites are divided in two classes according to whether the 
'local field' hi = —J'Y^^cos{6i — 6j — Aij) in the rt-th 

configuration 0^"^is greater or smaller than a threshold 
value ht chosen as explained below. The sites in the first 
group constitute the domains, ii) Correlations between 
the domains and the rest of the system are destroyed by 
a random rigid rotation of the former, iii) A fraction p 
of the iVw sites in weak local fields are picked at ran- 
dom and their phases reset to arbitrary values, iv) The 
energy of the subsystem formed by the domains is min- 
imized with the phases on the remaining sites fixed, v) 
The state resulting from the previous step is allowed to 
relax by performing a CGM of the total energy of the 
system. The outcome is the next state in the sequence, 
{6'n+i}. vi) The energy of this state is stored and, even- 
tually, hi is rescaled. 

The efficiency of this algorithm depends upon the 



chosice of the parameters ht and p. The threshold field 
fixes the degree of homogeneity required of a region for 
it to be classified as a domain. If it is too high or if p is 
too large, too many sites are involved in step iii) and the 
domain structure is disrupted just as in the SQA where 
all the phases are randomly reset at each step. 
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FIG. 1. Upper panel : Energies of some of the metastable 
states of a realization of the gauge glass model on a 512-site 
3D lattice. Results are plotted as a function of iteration num- 
ber for the spin quench algorithm (crosses) and for our algo- 
rithm (squares). The arrow indicates the energy of the first 
state in our sequence. Lower panel : Distribution of the en- 
ergies of two series of 5 000 metastable states each obtained 
with the SQA (a) and with our algorithm (b). 

If, on the contrary, ht or p are too small, the algorithm 
gets trapped in phase space and all the states in the se- 
quence are close to the initial one. The two parameters 
must therefore be continuously readjusted in the course 
of the simulation to ensure good performances. We have 
empirically found that the algorithm performs at its best 
for large samples when the number of sites involved in 
step in) pN^ = n ^ O.OSA^s- If at some stage of the it- 
eration < 71 we consider that the threshold field is 
too low and too many sites are being included in the do- 
mains. We then rescale it upwards, ht (l + a) ht with 
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a ~ 0.05, and we randomly reset the phases on all the 
sites where the local field is weak. If n < iVw < 2n, ht 
is unchanged and the phases are updated on just n ran- 
domly chosen sites. Finally, if Ny^ > 2n we consider that 
the threshold field is too high and we rescale it down- 
wards according to ht (1 — a) ht- Wc find that, in 
practice, the field stabilizes itself after a few iterations 
and oscillates about a value that, in the case of the sim- 
ulations reported here, is ^ 2.6J. 



In the upper panel of Fig, 



ll we show the energies of two 



series of 1000 minima of (Q) for a particular realization 
of the disorder. Data were obtained for a 3D lattice of 
512 sites using periodic boundary conditions. The crosses 
represent states obtained with the SQA and the squares 
are the outcome of the first thousand iterations of our al- 
gorithm. The arrow points at the energy of the first con- 
figuration of our sequence. It can be seen that, whereas 
the conventional algorithm randomly samples the whole 
of phase space, our method seems to mostly explore the 
deepest valleys. Notice that during the first five hun- 
dred or so iterations the typical energy of the states in 
the sequence decreases continuously after which it sta- 
bilizes in a region of energies that is hardly ever visited 
by the SQA. It is important to check that the configura- 
tions that enter in the sequence come from well separated 
regions of phase space rather than from a particular val- 
ley where the algorithm would be trapped. This may be 
done simply by monitoring the evolution of the overlap of 
the successive configurations with a particular one that 
is chosen as reference. The lower panel of Fig.|] shows 
histograms obtained after five thousand iterations of the 
two algorithms. It can be seen that the histogram ob- 
tained with our method is much narrower and centered 
at a much lower energy. The overall features of the dis- 
tribution of energies shown in Fig.|l| are quite simpar to 
those recently found for the ±J spin glass modelll. It 
is remarkable that about twenty percent of the states 
found using our algorithm in this example have never 
been generated by the SQA. Our lowest energy state, 
a.t E = — 2.036J, appears ^ 200 times in the sequence. 
The configurations of the states that have this energy are 
related to each other by uniform rotations. In between 
them, the algorithm generates states that are in far away 
regions of phase space. We believe that this state is the 
ground state of this particular realization. 

In order to study the scaling properties of defects ener- 
gies in the gauge glass model we have applied the above 
method to compute ground-state energies with periodic 
and antipcriodic boundary conditions for systems of 
sites with 3 < i < 8. OnlypSystems with L < 5 
had been investigated previouslySM . We generated se- 
quences of states containing 500 {L=3), 800 {L=A), 1000 
(L=5), 2000 (L=6), 3000 {L=7) and 5000 (i=8) ele- 
ments, respectively. Disorder averages were taken over 
25600 (L=3), 6400 (L=4), 2560 (£=5,6), 640 (L=7) and 
256 {L—8) samples, respectively. The normalized distri- 
butions of the differences AE — \Ep — Eap\ obtained 
numerically for the different sizes are shown in FigJ2. 



Detailed examination of the results shows that the dif- 
ferences between the curves for different sizes are of the 
same order of magnitude as the statistical error bars. Be- 
cause of this we were not able to determine the stiffness 
exponent of the model by performing a scaling plot as 
Eq.0 suggests. 
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FIG. 2. Probability distributions of the differences of the 
ground-state energies with periodic and antiperiodic bound- 
ary conditions. The solid line is the result of a gaussian fit of 
the ensemble of the data. 

The solid curve in the Fig.|^ is a gaussian of width 
W = 2.25J. The fact that we can quite reasonably 
describe the ensemble of the data using a single size- 
independent distribution is an indication that the lower 
critical dimension of the gauge glass pr|Qhlcj|a is very 

iMM. TheL- 
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close to three, as found by other authors 
dependence of the effective coupling W{L) 
is shown in the log-log plot of Fig.^. As the figure shows, 
statistics for the two largest systems is still unsatisfac- 
tory but very hard to improve upon because of CPU- 
time limitations. Nevertheless, we can still conclude from 
the available data that the domain- wall energy increases 
slowly with length scale. Leaving out the point for L — 3 
which is likely to be too small a size for scaling to hold, 
we can make a power-law fit of the results. The stiff- 
ness exponent thus determined is 9qq ~ 0.077 ±0.011. 
We have checked this result by repeating the calculation 
for the larger sizes starting from different random ini- 
tial configurations. The differences found between the 
results thus obtained fall within the statistical error bar. 
Our value for the stiffness, exponent is consistent with 
those reported Jav Gingraa3(6'GG — 0.04 ± 0.06) and by 
other authorsQ't3 who find 6qq « within their statis- 
tics. Ours is, to our knowledge, the first calculation in 
which the possibility Ogg < is outside the range cov- 
ered by the error bar. 



3 



J. R. Banavar and M. Cieplak, Phys. Rev. Lett. 48, 832 
(1982); M. Cieplak and J. R. Banavar, Phys. Rev. B 29, 
469 (1984). 

W. L. McMillan, Phys. Rev. B 31, 342 (1985). 

B. M. Morris, S. G. Colborne, M. A. Moore, A. J. Bray 

and J. Canisius, J. Phys. C 19, 1157 (1986). 

L. R. Walker and R. E. Walstedt, Phys. Rev. B 22, 3816 

(1980). 

P. Gawiec and D. R. Grempel, Phys. Rev. B 44, 2613 
(1991). 



1.0 1.2 1.4 1.6 1.8 2.0 2.2 
In/. 

FIG. 3. The L-dependence of the domain-wall energies for 
the gauge glass model on Lx Lx L simple cubic lattices. The 
dashed line is the power-law fit discussed in the text. 

The results of this paper thus confort the idea that 
the lower critical dimension of the gauge glass model is 
slightly below three. The implication is that the system 
has a finite-temperature transition to an ordered state in 
agreement wiih. the findings of the Monte Carlo studies 
of the modelBij. It is interesting to notice that whereas 
the smallness of ^gg would lead one to naively expect a 
very low transition temperature, the Monte Carlo data 
indicate that Tc ~ 0{J). This is a somewhat puzzling 
result that deserves further investigation. 

The calculations presented here have been done on a 
256-processor CRAY T3E parallel computer at the 'Cen- 
tre Grenoblois de Calcul Vectoriel'. We thank the staff 
for their technical help. 
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